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ABSTRACT 

This  study  is  to  understand  the  effects  of  nonlinear 
behavior  of  structures  using  a  finite  element  method. 

The  nonlinear  behavior  equations  are  derived  from 
equations  of  motion  and  constitutive  equations.  The  basic 
theory  is  the  principle  of  virtual  work. 

The  thesis  begins  with  a  comprehensive  formulation  of 
continuum-based  finite  element  theory.  The  theoretical 
portion  concludes  with  the  details  of  both  spatial  and 
temporal  discretization,  including  a  discussion  of 
nonlinearity . 

In  particular,  large  displacement  problems  and  visco- 
elastic  problems  remain  a  challenging  engineering  problems 
today.  The  viscoelastic  problems  depend  on  the  relaxation 
function  which  is  the  source  of  material  nonlinearities . 
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I.  INTRODUCTION 

A.  ANALYSIS  OF  STRUCTURES 

Analysis  of  structures  subjected  to  large  displacements 
requires  the  inclusion  of  nonlinear  terms  in  the  strain 
displacement  relations.  In  the  formulation  of  the  problem, 
either  total  Lagrangian  or  Eulerian  approach  is  adopted.  In 
order  to  assure  structural  integrity  for  loads  beyond  linear 
range,  nonlinear  analysis  is  very  much  mandated. 

Another  class  of  problems  is  the  viscoelastic  behavior 
of  structures  -  problems  where  stresses  in  the  components 
vary  with  time  at  a  constant  deformation  state  or  strains 
vary  with  time  at  constant  stress  states. 

These  problems  as  applied  to  simple  one-dimensional 
structures  will  be  addressed  in  this  research  work.  In  what 
follows  is  a  brief  survey  of  literature  in  this  area.  At 
the  end  of  this  chapter,  overview  of  the  thesis  is 
described. 

B.  LITERATURE  REVIEW 

A  method  for  nonlinear  static  and  dynamic  analysis  using 
cable  and  beam   elements   is   given   by   Bergan   et   al  [1]. 
Formulation   allows   for   material   nonlinearities   and  very 
large   deformations.     Gadala   and   Oravas   [2]   present   a 
comparison   of   finite   element   formulations   based   on  the 


Lagrangian,  updated  Lagrangian  and  Eulerian  formulations  for 
nonlinear  mechanics  problems. 

Henriksen  [3]  applies  a  single  integral  constitutive  law 
to  describe  viscoelastic  behavior  in  his  formulation. 
Ghoneim  et  al  [4]  treat  the  total  strain  rate  as  sum  of 
viscoelastic  and  viscoplastic  contributions  in  implementing 
the  constitutive  equations  in  their  formulation,  Pinsky  and 
Kim  [5]  use  variational  approach  for  studying  the  nonlinear 
viscoelastic  shell  behavior  by  representing  the  layered 
shell  as  a  multi-director  field.  The  formulation  includes, 
within  individual  layers,  the  effects  of  transverse  shear 
and  transverse  normal  strain  to  arbitrary  orders  in  the 
layer  thickness  coordinate  [6,  7]. 

Saigel  and  Yang  present  a  quadrilateral  shell  element 
and  describe  implementation  of  elastic-viscoplatic 
constitutive  law  for  large  deformation.  A  large  collection 
of  papers  on  the  application  of  ADINA  may  be  found  in  [8]. 

Response  of  nonlinear  mechanical  systems  to  random 
excitation  is  given  by  To  [9].  Thermomechanical  response  of 
uniaxial  bars  with  thermoviscoplastic  constitutive  laws  is 
discussed  by  Allen  [10].  Large  deformations  of  the 
contacting  bodies  on  the  contact  surfaces  and  solution 
methods  are  given  consideration  in  [11,  12].  Large 
deformations  of  composite  beams  subjected  to  axial  dynamic 
load  is  reported  in  [13]. 

Based   on   Eulerian-Lagrangian   formulation,   frictional 


contact  analysis  is  discussed  by  Haber  et  al  [14]. 

Halleux  and  Casadei  [6]  present  transient  large  strain 
analysis  of  structures  based  on  biquadratic  elements. 
Polar-decomposition  is  used  in  arriving  at  explicitly 
derived  tangent  stiffness  matrices  for  a  frame  element  by 
Kondoh  and  Atluri  [15].  A  method  based  on  co-rotational 
coordinate  system,  where  the  local  cartesian  coordinate 
system  attached  to  each  element  translate  and  rotate  is 
described  by  Mattiasson  [16]. 

Inelastic  problems  based  on  incremental  virtual  work 
equations  is  given  by  Nagtegaal  [17].  A  constant  strain 
rate  is  assumed  during  an  increment  in  integrating  the 
constitutive  equations. 

Planar  curved  beams  undergoing  large  rotations  is 
analyzed  using  Helinger-Reissner  mixed  variational  principle 
by  Noor  [18].  Steele  discusses  several  nonlinearities 
including  buckling,  creep,  and  nonlinear  elastic  problems  in 
[19].  The  use  of  Lagrangian  formulation  using  an  elasto- 
plastic  model  is  presented  by  Voyiadjis  and  Navaee  [20]. 
Hirashima  et  al  [21]  present  an  approach  for  the  snap- 
through  analysis  of  curved  shell  panels  and  determination  of 
initial  stress  matrix. 

Importance  of  viscoelastic  relaxation  effects  on  a 
pellet  are  investigated  by  Basombrio  [22],  where  it  is 
observed  that  the  duration  of  the  ramp  input  influences  the 
response  considerably.   Viscoplastic  constitutive  models  are 


obtained  explicitly,  without  resorting  to  inversion,  to  be 
directly  implemented  in  a  finite  element  formulation  [23]. 
An  intrinsic  time  scale  that  can  vary  with  time  and  position 
is  used,  together  with  means  to  control  the  numerical 
integration  in  the  approximation  of  constitutive  model  by 
Chambers  and  Becker  [24]. 

Chandra  et  al  [25,  26]  discuss  applications  of  visco- 
plastic  models  in  the  metal-forming.  Elasto-plastic  models 
are  used  in  analyzing  unilateral  contact  and  friction 
problems  by  Cheng  et  al  [12].  Chung  et  al  [27]  suggest 
methods  for  varying  parameters  in  guiding  solutions  for  the 
simulation  of  plastic  forming  processes.  Liu  et  al  [28] 
propose  formulation  for  path-dependent  materials  based  on 
arbitrary  Lagrangian  Euler  method  and  include  an  elastic- 
plastic  wave   propagation  problem. 

Rheinboldt  and  Riks  describe  several  algorithms  based  on 
continuation  method  and  mesh  refinements  for  the  solution  of 
nonlinear  finite  element  equations.  Ryu  and  Arora  [29] 
propose  an  approach  to  solve  nonlinear  finite  element 
equations  based  on  substructures  concept. 

Several  papers  pertaining  to  visco-elastic ,  visco- 
plastic  behavior  and  their  implementation  is  reported  in 
[10,  7,  4,  30,  31]. 


10 


C.   THESIS  OVERVIEW 

The  objective  of  this  thesis  is  to  review  some  of  the 
literature  available  in  the  analysis  of  large  deformation  of 
structures,  to  develop  the  theory  to  include  large 
displacements  and  viscoelastic  behavior  of  one  dimensional 
structures,  and  in  the  process,  to  understand  the  physics  of 
the  effects  of  nonlinearities .  Chapter  II  presents 
theoretical  formulation  based  on  principle  of  virtual  work 
and  incremental  load  method.  Descretization  based  on  a 
simplex  element  is  given. 

Chapter  III  gives  the  algorithms  used  in  the  solution 
procedure.  Full  Newton-Raphson  and  Modified  Newton-Raphson 
methods  are  described. 

Chapter  IV  describes  some  test  cases  and  numerical 
results  obtained. 

Finally,  Chapter  V  provides  conclusions  and  some  remarks 
for  future  work  in  this  area. 
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II.  THEORETICAL  FORMULATION 

A.  OVERVIEW 

In  this  chapter  we  consider  a  nonlinear  continuum 
equations  for  a  rod  element.  This  approach  is  based  on  the 
principle  of  virtual  work  [48].  The  resulting  equations  are 
in  convenient  form  for  numerical  solution  via  finite  element 
methods  [49].  This,  of  course,  is  the  conceptual  basis  of 
the  discretized-solid  finite  elements,  introduced  by 
J . T . Oden  et  al  [50]  in  the  elastic  structures  and  nonlinear 
continua. 

B.  NONLINEAR  CONTINUUM  EQUATIONS 

Discuss  a  rod  element  with  the  total  Lagrangian  Method. 
It  may  be  proved  that  internal  virtual  work  is  equal  to 
external  virtual  work.  Namely  a  deformable  system  is  in 
equilibrium  if  and  only  if  the  total  external  virtual  work 
is  equal  to  the  total  internal  virtual  work  for  every 
virtual  displacement  consistent  with  the  constraints. 

1.   Principle  of  Virtual  Work 

In  this  section  we  prove  that  total  internal  virtual 
work   is  equal  to  total  external  virtual  work.   Consider  the 
concept  of   minimum  potential  energy  [51].   The  condition 
is   satisfied   to   equilibrium   state   and    the   change   of 
potential  with   respect  to  the  displacement 
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must   remain   stationary.     And  consider  the  load- 

displacement  and  stress-strain  relationship        of  figures 
shown  below. 


-^suk- 


Figure  2.1  Load-Displacement  Curve 

6- 


Figure  2.2  Stress-Strain  Curve 


In   General,   the   total   potential   energy   is  equal  to  the 

summation  of  strain  energy  and  potential  energy. 

ie,    7TP   =   U   +   W  (2.1) 


where,    Kp 
0 

w 


Total  potential  energy 

Strain  energy  (internal  work) 

Potential  energy  of  load  system  (external 

work) 
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From  Figure  (2.1), 

W   =  -  J  p  6u 


=  -  2  Pi  ui 

where,    Pi  ith    load 

ui  deflection  under  ith  load 

From  Figure  (2.2), 


(2.2) 


U  = 


\  (\ 


o    6€  d(vol) 


vol   o 


5  J  a   €  d(vol) 

vol 


(2.3) 


Substituting  the  equation  (2.2)  and  (2.3)  in  equation  (2.1), 
we  get 


p   vol 


5  a  €  d(vol)  -  2  p.  u. 


(2.4) 


Equation   (2.4)   can   be   represented  by   matrix   form,  the 
following  is  the  matrix  formula. 


tup  =vJx-5  {a}    {€}  d(vol)  -  {  u.  }  {  p.  } 


2.5 


where,    {a} 
{€} 
(ui  } 
{Pi} 


a  vector  of  stresses 

a  corresponding  vector  of  strains 

a  vector  of  elemental  nodal  displacements 

a  vector  of  elemental  nodel  forces 


Let  us  define  [B] ,  a  matrix  relating  strains  and 
displacements , 


le 


{€}  =  [B]  {u} 


(2.6) 


and  define  [D] ,  a  matrix  relating  stresses  and  strains, 


le, 


{a}  =  [D]  {€} 


(2.7) 
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Substituting  equation   (2.6)  &   (2.7)  in   equation  (2.5),  we 
get 


Kp     - 


T      T  T 

{ui  }   {B}   [D]  [B]  d(vol)  {ui  }  -  {ui  }   {pi  }  (2.8) 


vol 

Now  we  consider  the  minimum  potential   energy,  the  condition 
is  as  follows: 

dTTp 


=  0 


(2.9 


d{ui  } 
Taking  derivatives   with  respect   to  {ui  }   in  equation  (2.8) 
and  applying  to  equation  (2.9),  we  get 


{pi  }  :J  [B]   CD]  [B]  d(vol)  {ui  } 
vol 


(2.10) 


From   relating   loads   and   displacements,   we   can   get  the 
elemental  stiffness  matrix  [k]  . 


le, 


[K]  =  J 


[B]   [D]  [B]  d(vol) 


(2.  11 


vol 


Summing  all   the  elements   of  the  entire  structure,  equation 

(2.10)  becomes 

{P}  =  [K]  {u}  (2.12) 

where,  {P}  -   1    {Pi } 
[K]   =  2    [K] 

[u]   =  2    {ui  } 


2 .   One  Dimensional  Finite  Element 

Consider  a  shape  function  which  is  linearly  function 
with   geometric   boundary   condition   with   one   dimensional 
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finite  element.   We  define  that  the  shape  function,   tn,  is 

\J/M(x)  =  aN   +   bN  X 

In  order  to  derive  the  constants  in  the  shape  function  for  a 
rod   element   with   linear    displacement   field,    use   the 
geometric  boundary  condition. 
We  may  obtain  the  constants  as  follows: 

X2 

i      Lo 


b    = (2.13) 

1  Ln 


XI 


a 


2  Lo 

1 


b 

2 


The  relationship  between  displacement  and  shape  function  is 
given  below 

M  T 

u   =  N|\{/N  uN  (N=l;  M=2,  ....) 

=  ^(x)  ux(t)  +  ^.(x)  u2(t)  (2.14) 

Another  consideration  is  the  viscoelastic  materials.  For 
viscoelastic  materials,  the  stresses  may  be  expressed  in 
terms  of  relaxation  functions  as  given  by 

't 


s(t)=  J      G(t  -  £  )  s    (£)  d£ 
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=  S   (t) 


I(t) 


(2.15) 


e  e 

where,    S   (t)  =  S   [ r  ( t ) ] 


I(t)  = 


G  (  t  -  £  ) 


ae 


:  elastic  response 
S*  (  £  )  d£ 


G(t):  relaxation  function. 
Assume  the   element  is  in  equilibrium  and  quasistatic  state. 
The  motion  equation  for  a  one  dimensional  finite   element  is 
given  by  J . T . ODEN  [ref.6], 


S\^     (  1  *£/  u   )  dV   =  P   (u  ,t 


Vo 


N  ,  x 


M,  x    M 


N     M 


where,     Vo   =  Ao Lo  :  undeformed  volume 

P    :  total  generalized  force  at  node  N 

N 


On   substituting    for   above    shape   functions   and   their 
derivatives,  we  may  obtain 

fN  (uM,  t,  I)  =  VobN  [  Se  -  I]  (1  +  bk  uk  )  -  PN(uM,t) 


=  0 


(2.16) 


C.   CONSISTENT  NODAL  FORCES 

Consider  a  rod  element  in  two  configurations:  undeformed 
and  deformed. 

Let  the  distributed  forces  per  unit  length  in  these 
configurations  be  q  and  q*  respectively.  The  consistent 
load  vector  depends  on  both  the  magnitude  of  the  load  and 
the  deformed  state. 
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One  dimensional  rod  element   under  distributed   load  is 
showed  as  following  figures. 


(x,V> 


-*. *. * 


i 


2- 


L. 


%*C*>*^ 


■■  i 


V 


a 


2,-- 


where,  Lo :  undeformed  rod  element  length 
L  :  deformed  rod  element  length 

q(x,t)  is  distributed  load  per  unit  undeformed  length 
q*(y,t)  is  distributed  load  per  unit  deformed  length 

Figure   2.3   One-Dimensional   Rod  Element  under  Distributed 
Load . 

Consider   the    relationship    between    the    deformed   and 
undeformed  coordinate.   It  may  be  written  as 

y(x,t)-x+u(x,t)  (2.17) 


where  x:  undeformed  coordinate  of  a  point 
y:  deformed  coordinate  of  a  point 
u:  displacement  of  a  point  within  element 

The  consistent  nodal  forces  vector  is  given  by 

In  index  nutation,  nodal  forces  may.be  written  as 


P  =  [  Pp%dV   -   fjh*  dVo 


(2.18) 


(2.19) 


with   N  =  1,  2 

where   F   :   body  forces  per  unit  mass  deformed 

Fo  :   body  forces  per  unit  mass  undeformed 

P     :   deformed  density 

P      :   undeformed  density 
o 

Relating  q*(y,t)  and  q(x,t)  is  given  by 

q(x,t)  =  \  q*(y,t)  (2.20) 

where,   A.  =  dy/dx     :  stretch 

Consider  the   undeformed  body   force  per   unit  mass,  then  we 
may  write 


r°"    dm     -     P0Ao     "      P.  A 


(2.21) 


Using  equation  (2.17) 

y   =  x  +  \Jf   u  (2.22) 

M   M 

dy 
=  1  +  \|/     u  (2.23) 

dx  M,  x    M 

The  undeformed  body  force  per  unit  mass  is  then 

A      q* 

P  =  (  1  +^M,x  um)  (2.24) 

'  0  o  Ao 

The  consistent  nodal  forces  take  the  following  form, 

fL     t       fL° 
PN  =  J0  q*  Vf?x  +  Joq*>^MiXu  M\^Ndx  (2.25) 

Substituting  the  shape  function  in  equation  (2.25),  we  may 
get 
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PLo 

5N  =  (1  +  bMuM)  J    q*(aN  +  bN  x)  dx 


(2.26) 


Now  we  consider  some  special  cases  of  equation  (2.26) 
One  case  is  the  uniformly  distributed  loads. 


/ 


S* 


V 


*>* 


Figure  2.4  Rod  Element  under  Uniformly  Distributed  Loads 


The  uniformly  distributed  loads  are  given 

q*(y,t)  =  qo   h(t) 
Substituting  the  loads  in   equation  (2.26), 
consistent  nodal  forces. 


we  may   get  the 


{  P„  >  = 


where,    h(t 


p  /I    qo  h(t) 

>  = [Lo  + 


X    J 


u*-ui]U 


du 


dx 


Another  case  is  the  linearly  distributed  loads 


(2.27) 


Figure  2.5  Rod  Element  under  Linearly  Distributed  Loads 
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Linearly  distributed  loads  are  given 


q(x,t)  =  q  (t)  +  [q  (t)  -  q  (t)] 


ey 


Substituting   the   loads   in   equation  (2.26)  we  may  get  the 
consistent  nodal  forces. 


<  PN   > 


13 


(Lo  -  u   +  ua  ) 


l" 

M 

2 

< 

,<\ 

»   (2.28) 


D.   SOLUTION  METHODS 

In  this  section  we  derive  the  equilibrium  equations  of 
incremental  method,  Full  Newton-Raphson,  and  Modified 
Newton-Raphson  Method. 

It  is  based  on  total  Lagrangian  Method  and  use  the 
Taylor  Serier  at  a  given  equilibrium  point. 

The  element  stiffness  matrices  and  consistent  load 
vectors  are  assembled  to  form  global  matrices  in  the  usual 
manner.  The  prescribed  displacements  are  implemented  by 
replacing  the  diagonal  element  of  the  corresponding  row  by 
unity  and  the  rest  of  the  element  of  the  rows  and  columns  by 
zeros  [49] . 

1 .   Incremental  Method  (INC) 

Consider  the  incremental  process,  we  can  rewrite  the 
equation  of  load-displacement  relationship, 

[K]i  (A>i+i  r  (APMi  +  {APc}i  (2.29) 

where,      i  :  the  current  configuration  coordinates. 
{ APc  }i  :  the  force  unbalance 
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{/\P}i+i  :  an  increment  in  the  external  load 
The  method  is  described  by  the  Figure  2.6 


Pb 

p. 


c 

c 

-A 

,  ....v*    ' 

\ 

H 

/f    "ICj 

1 

1 

1 

A 

JL 

As         *s>  1         J 

'          \                 1 

-p. 

1 

1 

k    «■  * — ^—4 

1 

1 

where,  AEC    The  algorithm  of  equation  (2.29) 

AB'  C  :  The  algorithm  of  purely  incremental 
process 

Figure  2.6.  The  Algorithm  of  Incremental  process 


The  approximate  values  of  A,  E,  C,  is  computed  from  p  -f(u). 
To  improve  these  values,  take  the  iteration  process  at  each 
load  level  with  {AP}i+i=  0,  as  is  Figure  2.7. 

The  iteration  processes  are  two  ways,  one  is  few 
load  levels  with  many  corrective  iterations  in  each,  or 
another  is  many  load  levels  with  few  iterations  in  each. 
This  is  sometimes  called  'incremental  with  one-step  Newton- 
Raphson  method" 

2.   Full  Newton-Raphsor:  Method  CFNR) 

Consider  the  load  versus   displacement  function   P  - 
f(u).   Take  the  Taylor  series  of  P  =  f(u) 
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JCUA+   A,)=  J(U^  ♦  <#>.    A, 


(2.30) 


dP 


where , 


K 


tangent  stiffness  at  A 


du 


f(u  ) 

A 


We 


f  (u   +  A  )  =  P 


look  for  A.such  that   f(u   +  A  ) 


-     P   ,   SO 
B 


k     (A)  =  p    -   P 

A       I       B       B 

The  consequent  steps  are  explained  by  the  Figure  2.7 


(2.31) 


Figure  2.7.  The  Newton  Raphson  solution  of  the  equation 
P  =   f(u) 

The  above  algorithms  are  below; 

update  displacements  :  ui  =  l|A  +  A , 

use  ui  to  get  stiffness  ki  and  load  Pi 
.  find  next  increment  ^a  from   ki  (  Ax)  =  Pb  -  Pi 

Eventually   ua   +   A  ,    +•   A  ^      +  =  ub  to  a  close 

approximation . 

Then,  the  load  can  be   increased   from   Pb   to   Pc ,   and  the 

process  can   be  started   again.    This  solution  algorithm  is 
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known  as  the  Newton-Raphson  Method. 

3.   Modified  Newton-Raphson  Method  (MNR) 

The  Newton-Raphson  method  requires  that  tangent 
stiffness  matrix  [k]  be  generated  and  reduced  for  equation 
solving  in  every  iteration  cycle.  This  process  is  not 
economic.  An  alternative  is  modified  Newton-Raphson 
iteration,  also  called  "constant  stiffness  iteration".  Here 
the  same  tangent  stiffness  matrix  is  used  for  several 
iterations  cycles.   The  MNR  is  described  by  the  Figure  2.8 


Figure  2.8  The  MNR  Iteration 

As  compared  with  FNR,  the  MNR  requires  more  iterations, 
but  each  is  done  more  quickly  because  [K]  is  formed  and 
reduced  in  only  the  first  iteration. 

E.   ANALYTIC  SOLUTIONS 

At  first  we  define  the  concept  of  stresses  [52]  Eulerian 
stress,  a    ,    is  defined  as  follows: 
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N(y,t) 

a    (y,t)  =  (2.32) 

A(y) 

and  is  called  True  stress.   Lagrangian  stress,  T,  is  defined 

as  follows: 

N(y,t) 

T(x,t)  -  (2.33) 

Ao  (x) 

and  is  called  Engineering  stress. 

Kirchloff  stress,  S,  is  defined  as  follows: 

N(y,t) 

S(x,t)  =  — (2.34) 

\(x,t)  Ao(x) 

The  relationship  between  strain   r  ,  and   stretch      is  as 

follows : 


(  \*  "  1 ) 


du 

\  =  i  +  (2.36) 

dx 

where,    N  :  Applied  force 

A  :  deformed  area 


For  the   special  case   of  small  strain,  both  u,  x  and  change 

in  cross-sectional   area  are   small.    In  this   case,  we  can 

write 

r  *  u   =  € 
a   *  T  »  S 

where,  €   is  engineering  strain. 

Thus  in   small  deformation  case,  the  Lagrangian  stress  T 
and  engineering   strain   €     can   be   used   to   replace  the 
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and   engineering   strain   €     can   be   used   to  replace  the 
Kirchhoff  stress  S  and  Lagrangian  strain   r  . 
1.   Elastic  Materials 

Consider   a   one-dimensional   rod   with  varying 
cross-  sectional   area,  Ao  =  Ao(x),  in  equilibrium  under 
distributed  external  loading. 


AoOO 


(T+£*X) 


(A*+%  *x) 


(a)  (b) 

where,   f    external  distributed  load  per  unit  length 

T  :  Lagrangian  stress  within  rod 
Figure  2.9.  One-Dimensional  Rod  under  Axial  Loading 


For  static  problems  and  elastic  materials,  both  f   and  T  are 
independent  of  time, 


£       —       -C  f  ..   \       n«-J  1*       -       rr'  i   •  -  x 


From  figure   2.9(b),  the  static  equilibrium  of  forces  in  the 
x-direction  is 

(2. 37 ) 
£  F(x)  =  -TA0  +  fax  ♦  (T  ♦  (dT/dx)  ax  )  (A0  ♦  (dA0  /  dx)  Ax  )  =  0 

Expanding  in  Taylor  series, 
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^c  =  0 


(2.38) 


For  a  continuous  rod, 


in  the  limit, 


i+i^TA^  - 


(2.39) 


Substitution  of  equations  (2.34),  (2.35),  and  (2.36)  in  the 
above  (2.40) 

where,    u(x)  :  displacement  w.r.t.  x 

u, x   :  first  derivative  of  u(x) 

u.xx  :  second  derivative  of  u(x) 
Note  that  equation  (2.40)  gives  the  values  of  distributed 
force  f  (per  unit  length)  required  to  produce  the 
displacement  field  u(x)  of  a  one-dimensional  rod  once  the 
boundary  conditions,  the  stress-strain  law,  and  varying  area 
function  Ao ( x )  have  been  specified.  Consider  the  rod  with 
constant  area  Ao ,  equation  (2.40)  becomes 


f(x)  =  -A0  u,xx  1(1  +  u,x )  **2  (  dS*  /dr )  +  S«] 


(2.41) 


Some  special  cases  of   equation  (2.41)   is  considered.    One 
case  is  the  specified  displacements. 
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/- 

/- 

/- 


\oo> 


-=*  X 


JA« 


where,  u(x)  -    ao  +  a, x  +  aa x2 

u,  x  "  ai  +  2  aa  x 

u,  x  x  =  2  a» 
Figure  2.10  Test  Model  of  "Inverse  Solution" 

Apply  the  boundary  contition  u,x(x)  =0   at  x  =  Lo 

a 

u  (  x  )  =  a  ,  x  — -  x2  (2.42) 

2Lo 

If  at   is   specified,  we  may  get  the  exact  solution  u(x)  and 

the  corresponding  required  distributed  load  f(x). 

a.   Linear  displacements  field: 

u(x)   =   a, x 

u,  x   =   a, 

U,  xx   =   0 


From  equation  (2.42),  we  may  get 

f(x)   =  0 
Note  that  a  linear   displacement   field   cannot   be  produced 
with  distributed  load. 
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b.   Quadratic  displacement  field: 
u(x)   =  a,x  -  (ai/2Lo)x2 
u,  x   =  a{  -    y  (  a  i  /  Lo  )  x 
11,  u  -    -  ( a  |  /  Lo  ) 

-  Consider  the  linear  materials,  Se  -    E  r 

We  may  get  the  distributed  forces  from  equation  (2.41) 
f(x)  -      -    EAuU.xx   [(1  +  u,  x  )  +  (u,  x  +  %   u2,  x  )]     (2.43 

-  -  Consider  the  nonlinear  materials, 
Se  =  Ei r  +  E2 r2  (quadratic  form) 
From  equation  (2.41),  we  may  get 

f(x)  =  -  AoU.xx  {  (1  +  u,  x  )  [Ei  +  2  E2  (u,x  +  %   u2,x)] 

+  El  (u,  x  *-  ^u2,x)  +  E2  (u,x  +  >2U2,x)]  }     (2.44) 
Another  case  is  the  prescribed  loads. 


/ 


/ 


A 


e 


^p 


-*  3C 


Figure  2.11  Test  Model  of  Rod  with  Prescribed  End  Element 


From  equation  (2.39),  since  no  distributed  loads  exist, 

d 

( TAo )  =  0   with  constant  area. 

dx 

We  may  get 

G  =  constant  (2.45) 
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The  applied  load  is  given 

P  -  Ao\  S  (2.46) 

-  Linear  material; 

^e  r  Er  =  E(u,  x  +  *i  u^  ,  x  ) 
From  equation  (2.46),  we  may  get 

P  =  EAo  (u,  x  +    u2,x  +  ^  u3,x)  (2.47) 

-  Nonlinear  material; 

Se   =   El r   +   E2 r2   (quadratic  form) 
From  equation  (2.46)  we  may  get 

P  =  Ao   (1  +  U,  X  )   [El   (U,  x   +  y2    U=  ,x  )  +  E2   (U,  x   +  %    U2  ,x   )2] 

(2.48) 
2 .   Quasi linear  Viscoel a stic  materials 

The  Viscoelastic  materials  has  the  more  complicated 
problems  than  the  elastic  materials,  since  the  external 
loads  and  displacements  are  time-dependent. 


< 


lU-t) 


x^ 


u< 


uLoc*> 
>  Fu0c*) 


X=Lo 
— ^1 


where,   Ao ,  Lo  :  Underf ormed  Area  and  Length 

uo(t)  :  time-varying  left-end  displacement 

uLo(t)     time-varying  right-end  displacement 

Fo  ( t )  :  time-varying  left-end  force 

FLo(t)     time-varying  right-end  force 

Figure  2.12  One-Dimensional  Rod  Subjected  to   Time-varying 
End  Loads 
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For  quasilinear  viscoelastic  materials, 


Sit) 


Jn 


0  &, 


(2.15) 


The  Lagrangian  stress  is  constant  throughout  the  rod  since 
the  distributed  loads  is  not  applied.  From  equation  (2.39), 
the  stress  is  time  variable  only. 

T(x,t)  -  g(t)  =  \(t)   S(t)  (2.49) 

And  the  internal  forces  produced  with  the  rod  is 

F(t)   =  TAo  =   Aog(t)  (2.50) 

Substituting  equation  (2.15)  in  equation  (2.49)  yields 

g(t)   =  \(t)   S(t) 


(1  +  u.) 


Se(r)  - 


o  tf 


(2.51) 


The  right   hand  side   of  equation   (2.51)  is  a  time  function 

only . 

Therefore   g(t)   must  also  be  a  time  function. 


du 
dx 


u,x(x,t)  =  h(t) 


(2.52) 


Integrating  w.r.t.   x   at  equation  (2.52) 

u=   u(x,t)  -      xh(t)  +  C(t)  (2.53) 

where,  C  :  integrating  function  to  be  determined  by  the 

boundary  condition  at  one  end. 
Assume  that 

u  (  o  ,  t )  =  uo  ( t )  =  C  ( t )  , 
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therefore 

u(x,t)=xh(t)+uo(t),  (2.54) 

u,x (x, t)  =  h(t) 
Substituting  equation  (2.54)  in  equation  (2.51)  yields 


;] 


(2.55) 


Once  the  G(t)  and  Se ( r )  have  been  specified,  the  equation 
(2.55)  can  be  used  with  the  boundary  conditions  to  find  the 
exact  solution  in  some  cases.  Consider  the  3  element  solid 
with  linear  material  property.  The  prescribed  displacements 
boundary  conditions  at  both  ends  is  described  in  Figure 
2.13. 


/ 
/ 

/ 

/ 


/£ 


L, 


■>  x 


-J 


u(o, t)  =  uo ( t)  -  0 
u  (  Lo  ,  t )  =  u,  ( t ) 


(2.56) 


Figure  2.13  Test  Model  of  a  Rod  with  Prescribed  Time  Varying 
Displacement  at  the  End. 


The  prescribed  displacements  are  linear, 
u(x,  t)  -  (x/Lo  )  u.  (t) 


U,  X  (  X  ,  t  ) 


u.   ( t )  /  Lo  =  h  ( t ) 

I'  o 


(2.57) 
(2.58) 
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where,   t^  (t)  =  Loh(t) 

Assume  linear  material  property,  Se  =  Er 

And  for  a  3  element  rod,  G(t)is  defined  below 

G(t)  =  (1  -  a    )  +    e-Pt  (2.59) 

where,   a   ,  P  ;   constants  for  material  properties. 

Note  that  G(o)  =  1  at  t  =  0 ,   G(t)  =  1  -  a      at  (3-0 

No  relaxation  occurs  and  the  material  exhibits  elastic 
response . 

Specifying  the  prescribed  displacement  function  u  r  (t). 

L  o 

a.   Ramp  Function 

At  x  =  Lo  ,   B.  C.'S:u(o,t)=0 

u(Lo,t)=ru   (t)  =  at  (a:  constant) 

u  o 

From  equation  (2.57)  and  2.58),  we  may  get 

u(x,t)  =  (a/Lo )  x  t  (2.60) 

u,  x  (  x  ,  t )  =  (  a/Lo  )  t 
Substituting  equations  (2.57)  -  (2.59)  in  equation  (2.55) 
We  may  get  the  reaction  force,  F. 


EAo         at          a2 12    a    a       a     -|3t 
[1  +  ]  [{at+  } {!-( )}e 


Lo  Lo  2Lo      (3      pLo 

a  a  a 


{(Pt  -  1)  + (p2t2  -  2(3 1  +  2)}]    (2.61) 

(3  2  pLo 

Equation  (2.61)  describes  the  reaction  force  produced  at  the 

fixed  end  for  the  prescribed  displacement 

u  ,  ( t )  -    at 
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b.   Step  Function  (Relaxation  Test) 

At  x  =  Lo ,  B.C.'S  :  u ( o , t )  =  0 
u  ( Lo  ,  t )  =  ui^l  ( t ) 
where,  1 ( t )  :  unit  step  function 
From  equation  (2.57)  and  (2.58),  we  may  get 

ULo 

u(x,t)  =  x  l(t)  (2.62) 

Lo 

ULo 

u,  x i  (t)  (2.63) 

Lo 

Substituting  equations   (2.57)  -  (2.59)  and  (2.62)  in 

equation  (2.55),  we  may  get  the  reaction  force,  F. 


EAo  ULo  ULo 

F(t)  =  ( )  ulo  (1  +  )  G(t)  [1  +  ( )   1  (t)] 


(2.64) 


Lo  2Lo  Lo 

Consider  the  Kirchhoff  stresses. 
From  equation  (2.34) 

ULo  ULo 

S  =  E  ( )  (  1  +  )  G(t)  =  E  ru   G(t)      (2.65) 

Lo  2Lo 


ULo  ULo 

where,   rL   =  (  1  +  )     Lagrangian  strain  due  to 

Lo  2Lo 

a  prescribed   step  displacement  function  u,  l(t)  at  the  rod 

end. 
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III.  PROGRAM  IMPLEMENTATION 

A.   OVERVIEW 

We   are   interested   in   static   analysis.     For  static 
analysis,  load-displacement  relationship  is  as  following: 

(P)  =  [k]  {u} 
where  {P}    external  loads  vector 

[k]  :  constrained  stiffness  matrix 
{u}    nodal  displacements  vector 
If   we   know   {P}   or   {u},   then   we  can  analyze  the  static 
problems  by   using  elementary   elastic  theory.    In  general, 
the  static  analysis  steps  are  shown  as  follows: 
Step  1    Enter  geometry  of  material  properties,  loads, 

et  al . 
Step  2    Calculate  element  stiffness  matrices. 
Step  3    Determine  appropriate  part  of  [k]  and  {P}. 
Step  4  :  If  [k]  or  {P}  is  appropriate,  then  repeat  step  2 

through  step  3 . 
Step  5  :  If  [k]  or  {P}  is  not  appropriate,  then  apply 

constraints . 
Step  6  :  Determine  {u}. 
Step  7  :  Regenerate  elements. 
Step  8  :  Calculate  element  stresses. 
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B.   SOLUTION  PROCEDURES 

In  nonlinear  structural  problems,  to  find  the  solutions 
we  use  the  INC  and  FNR  and  MNR.  In  fact,  it  is  hard  to  know 
in  advance  what  increments  of  loads  or  iteration  should  be 
used  to  obtain  a  good  approximation  to  the  solution.  The 
interactive  method  was  used  in  determining  the  strategy  for 
solution  in  program  application.  The  strategy  subroutines 
were  developed  for  elastic  or  visco-elastic  time-dependent 
structural  problems.  Here  the  initial  state  for  incremental 
loading  process  was  selected  as  a  null  displacement  vector. 
Next  we  consider  finite  element  method  for  elastic  (FEMEL) 
and  viscoelastic  (FEMEV)  programs.  IN  both  programs, 
equilibrium  was  checked  at  each  stage  for  solutions  obtained 
from  other  INC  or  FNR  and  MNR. 

1 .   Algorithm  Implementation 

Consider   the   algorithm   of   Newton-Raphson   Method 
based  on    total   Lagrangian   Method.     Full  Newton-Raphson 
Algorithm  is    as  follows: 
Step  1.  Establish  local  coordinates  X-Y  for  the  element  at 

hand. 
Step  2.  Generate  the  element  stiffness  matrix  [k]  in  local 

coordinates  so  that  it  operates  on  the  local  D.O.F. 
Step  3.  Transform  [k]  to  global  coordinates  so  that  it 

operates  on  global  D.O.F. 
Step  4.  Repeat  step  1  through  3  until  all  elements  have  been 
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treated. 

Step  5.  Compute  displacements. 
Step  6.  Compute  internal  forces,  [k]{u}. 
Step  7.  Solve  [k]  {Au}  =  (Ap> 
where,  uN  +  1  =  uN+  A  uN 

Step  8.  Test  for  convergence. 

If  not  satisfied,  return  to  step  1. 

Where,  convergence  =  ((ARTAR)  /  (AprAp))**0.5  < 

0.001 
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2.   Flow  Chart  of  the  Finite  Element  Program 


F" 


<- 


In  put:  geometry, mat. prop;  loads 


3L 


PROBLEMTYPE: 
ELASTA/ISCOELASTIC 


ELE. MATRICES 
STIFF/LOAD: 


GLOBASSEMBLY 


SOLN:FORWARD 
ELIMINATION 


-r- 


SOLN:BACKWARD 
SUBSTITUTION 

£ 


UPDATE  DISPS 


4/ 


COMPUTE  RESIDUAL 
FORCES 


SOLUTION 
MONITOR 


DIAGNOSTICS 


NEW  LOAD  INC/UPDATE  DISPS  e 


PRINT  RESULTS 


STOP 


[Flowchart  of  the  Finite  element  Program  I 


IV.  NUMERICAL  EXAMPLES 

We   consider   two   numerical   examples   for   elastic  and 
viscoelastic  materials. 
A.   ELASTIC  MATERIAL 

1  •  Prescribed  concentrated forces . 

Specify  the  rod  element. 


/ 


/ 
/ 
/ 


L 


M 


»    Pz 


Pii  -  22,5000  lbs 
A   =  0.25  [in: j 
L   =10   [in] 
E  -    106  [lbs/in2] 

Consider  the  stiffness  matrix  of  rod  element,  [k]  . 
the   rod  element, 


For 


[K] 


EA   r  Kn    ki  2- 
L    |_k2i    kc  2 

Consider  the  boundary  conditions, 
ui   =  0 
Pi   =  0 


(4.1) 
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P2   =  22,5000  at  x  -  L 

From  loads  -  displacements  relationship, 

Solve  for  the  displacements,  u2  , 
Hence,  we  may  get 

LP2 

EA 


U2 


9.0  (in) 


2 .   Prescribed  displacements 


/ 


/ 


/ 


L 


-I 


X 


Consider  the  same  rod  element  before, 

the   quadratic  displacements,  u(x), 

u(x)   =   0.02  x   -   0.001  x  2 

The  first  derivative  is 

du 

—   =  0.0  2  -  0.002  x 

dx 

The  second  derivative  is 

d'^u 

=   0.0U2 

dx2 


and  specified 
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From  equation    (   2.41  )  with  the  constant  area,  we  may 
getthe  distributed  loads  for  any  x. 


d2 u         du     dse 

f  (  X  )   =  -  Ao  [  (  1  + )  2 +  S*  ] 

dx2         dx     dr 
=  3.0  *  10-3  X2   -   3.06  x   +   530  (4.2) 
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B.   VISCOELASTIC  MATERIAL 

Consider  the  rod  element  below, 


F.C*) 


U.(*> 


->  -ya 


Boundary  conditions;  u(0,t)  =  uo  (t)  =  0 

u(Lo,t)  =  u    (t)  =  u    1 ( t ) 
L  o  L  o 

Lagrangian  Stress;  T(x,t)  =  g(t)  -  A(t)  S(t) 
Internal  Force;  F(t)  =  T  Au  g(t) 


Kirchhoff  stress;    S  = 


(4.3) 


X(x  ,  t)  Ac  (x) 
For  3  element  solids,  the  relaxation  function  is 
G(t)  =  (1  -  a)  t  a  e-pt  (4.4) 

where,  OL  and   |3   are   constants   of   the   material 

properties . 

Assume  unit  step  function  for  prescribed  displacements. 

At  that  time,  the  boundary  condition  is 
u(o,  t)  =  0 
u  ( Lo  ,  t )  =  u.  ( t )  =  u.   1  ( t ) 

u  o         u  o 
where,  l(t):  unit  step  function 

We   may   get   the   displacements   function   with   unit   step 

function . 


u  (  x  ,  t )  =  -- -   ul  o  1  (  t 

Lo 


(4.5) 


4  5 


du     ULo 

—  =  l(t)  (4.6) 

dx     Lo 

And  the  Kirchhoff  stress  is  specified 

S   =  E  r  (4.7) 

Substituting  equations  (4.4)  -  (4.7)  in  equation  (4.3), 

EAo  UL  o  UL  o 

F ( t )  = Ulo  1  +  G  ( t )   1  + 1  ( t )   (4.8) 

Lo  2Lo  Lo 

Hence,  we  may  get  the  Kirchhoff  stress,  S, 

S  =  rr    G(t)  (4.9) 

L  o 

ULo  ULo 

where,   rLo  = 1  +  :  Lagrangian  strain 

Lo         2Lo 

due  to  the  prescribed  step  displacement  function 
u  1 ( t )  at  the  rod  end . 

Equation   (4.9)   describes   the  relaxation  cases.   From 
equation  (4.8),  we  may  get  the  internal  forces. 
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V.  CONCLUSION 

This  study  is  directed  towards  understanding  the 
nonlinear  behavior  of  structures. 

The  formulation  is  based  on  the  principle  of  virtual 
work.  We  consider  the  one  dimensional  nonlinear  continuum- 
based  rod  element  that  includes  large  displacements  and 
viscoelastic  effects.  To  get  the  solution  we  use  the 
incremental  process  with  full  Newton-Raphson  and  modified 
Newton-Raphson  Method.  Numerical  solutions  agree  well  with 
exact  solutions. 

Further  studies  may  be  extended  to  two  dimensional 
structures,  such  as  laminated  composite  plates. 
Consideration  of  the  viscoelastic  effects  of  composite 
panels  at  high  temperatures  is  another  possible  area  of 
study. 

An  experimental  study  to  verify  the  results  of  the 
present  study  is  also  recommended. 
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